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A set of key properties for an ideal dissipation scheme in gyrokinetic simulations is proposed, 
and implementation of a model collision operator satisfying these properties is described. This 
operator is based on the exact linearized test-particle collision operator, with approximations to the 
field-particle terms that preserve conservation laws and an H-Theorem. It includes energy diffusion, 
pitch-angle scattering, and finite Larmor effects corresponding to classical (real-space) diffusion. The 
numerical implementation in the continuum gyrokinetic code GS2 is fully implicit and guarantees 
exact satisfaction of conservation properties. Numerical results are presented showing that the 
correct physics is captured over the entire range of collisionalities, from the collisionless to the 
strongly collisional regimes, without recourse to artificial dissipation. 



I. INTRODUCTION 

Collisions play an important role in gyrokinetics. An 
accurate collision ope rato r is important for calculation of 
neoclassical transport!^ and the growth rate of instabil- 
ities s uch a s trapped electron modesPH dissipative drift 
wavesj^^ and microtearing modes^ in moderate colli- 
sionality regimes. Collisions can also affect the damp- 
ing of zonal flows^ and other modes that provide a sink 
for turbulent energy. In their a bsence, arbit rarily fine 
scales can develop in phase spacef^ 11 * 12 * 13 * 1 ^ which can 
in some cases pose challenges for discrete numer ical al- 
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gorithms, especially in the long-time limit 
modest amount of collisions can make accurate numeri- 
cal calculation much easier. 

Furthermore, inclusion of a small collisionality keeps 
the distribution function smooth enough in velocity space 
that the standard gyrokinetic ordering^ for velocity 
space gradi ents i s satisfied. For example, the parallel 
nonlinearityj 18 * 19 ! given by 
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enters at the same order as the other terms in the gy- 
rokinetic equation if the typical scale of parallel veloc- 
ity fluctuations, 8v\\, is one order smaller in the gy- 
rokinetic expansion parameter p/L (p = gyroradius and 
L = background scale length) than the thermal speed, 
Vth- Here, h is the non-Boltzmann part of the perturbed 
distribution function (defined more rigorously in the next 
section), ip is the electrostatic potential, B is the mag- 



netic field strength, b = B/_B, q is particle charge, in is 
particle mass, and ( . ) denotes the gyroaverage at fixed 
guiding center position R. 

While such a situation is possible in the collisionless 
limit, a small collisionality prohibits the formation of 
structures with 6v\\ ~ (p/L)v t f l . The level of collision- 
ality necessary to negate the importance of the parallel 
nonlinearity can be calculated by assuming a balance be- 
tween collisions and fluctation dynamics: 
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where C[h] describes the effect of collisions on h, u> is the 
fluctuation frequency, and v is the collision frequency. 
From the above expression, we see that scales in velocity 
space become small enough for the parallel nonlinearity 
to be important only when the collision frequency satis- 
fies v ~ (p/L) 2 uj. Such low collisionalities are not present 
in most fusion plasmas of interest. Furthermore, if such 
an ordering had to be adopted, the lowest order distri- 
bution function could become strongly non-Maxwcllian. 
This is a problem for Sf codes that assume an equilib- 
rium Maxwellian. 

In light of the above considerations, it is important to 
include an accurate treatment of dissipation in gyroki- 
netic simulations. In order to faithfully represent gyroki- 
netic plasma dynamics at reasonable numerical expense, 
we take the view that the form of the dissipation should 
be such that it: ensures satisfaction of the standard gy- 
rokinetic ordering; locally conserves particle number, mo- 
mentum, and energy; satisfies Boltzmann's _ff-Theorem; 
and efficiently smooths phase space structure. The first 
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of these requirements has already been discussed in the 
context of the parallel nonlinearity. Conservation prop- 
erties have been found to be important, for instance, in 
calculations of the neoclassical ion thermal conductiv- 
ity,^ as well as in a wide range of problems in fluid dy- 
namics. The e xistence of an ff-Theorem is critical for 
entropy balance^ * 21 * 22 ! and for th e dynamics of the tur- 
bulent phase space cascade? 12 * 13 ! Efficient smoothing of 
phase space structures is necessary to resolve numerical 
simulations at reasonable computational expense. 

A commonly employed dissipation mechanism in gy- 
rokinetic simulations is artificial (hyper) dis sipation, of- 
ten in physical (position) spaceP ^ | 23 | 24 | 2S | Ideally, the 
form of the artificial dissipation should be chosen to sat- 
isfy the requirements listed above and should be tested 
for convergence to the collisionlcss result. Of course, arti- 
ficial dissipation alone is unable to capture the correct dy- 
namics for moderate to strongly collisional systems where 
turbulent fluxes and other observable quantities depend 
sensitively on collisionality; for such systems, a physical 
dissipation model is desired. 

A number of such model physical collision operators 
are employed in gyrokinetic codesP^ 3 * 24 ! 26 ! These range in 
complexity from the Krook operator^ to the Rutherford- 
Kovrizhnikh operator^ to the Catto-Tsang operator,^ 
all of which have previously been implemented in GS2 
(see, e.g. Ref. 20). However, none of these satisfy all 
of the properties we require of a good collision opera- 
tor (See Appendix C of the accompanying paper^S! for 
a fuller discussion of this point). Here, we discuss nu- 
merical implementation in GS2 of an improved model 
operator which: includes the effects of both pitch- angle 
scattering and energy diffusion (i.e. efficiently smooths 
in phase space and ensures gyrokinetic ordering); con- 
serves particle number, momentum, and energy; satisfies 
Boltzmann's ff-Theorem; and reduces to the linearized 
Landau test-particle operator in the large k±p limit. A 
full description of this operator and a discussion of its 
desirable properties, including comparison with the pre- 
vious models of Catto-Tsang and Hirshman-SigmarpD is 
given in the companion paper 3 -^ (henceforth Paper I). 
We will focus on how such an operator can be imple- 
mented efficiently in gyrokinetic codes while maintaining 
the properties listed above and on how our gyrokinetic 
dissipation scheme (or any other) might be tested against 
a number of plasma physics problems. 

The paper is organized as follows: in Section II, we 
introduce the gyroaveraged collision operator from Pa- 
per I and examine properties that should be taken into 
account when using it in numerical simulations; in Sec- 
tion III, we describe our numerical implementation of the 
collision operator; in Section IV we present numerical re- 
sults for a number of tests demonstrating the ability of 
our collision operator implementation to reproduce cor- 
rect collisional and collisionless physics; and in Section 
V, we summarize our findings. 



II. PROPERTIES OF THE GYROAVERAGED 
COLLISION OPERATOR 

In order to include collisions in gyrokinetics, we fol- 
low the treatment of Ref. and assume the collision 
frequency, to be the same order in the gyrokinetic or- 
dering as the characteristic fluctuation frequency, wP^l 
This leads to the requirement that the distribution of 
particles in velocity space is Maxwellian to lowest order 
and allows us to represent the total distribution function 
through first order in p/L (where p is ion gyroradius and 
L is the scale length of equilibrium quantities) as 



f(T,n,e,t)=F (e)[l 
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where r is particle position, R = r bxv/f^o is guid- 
ing center position, p = mv\/2Ba is magnetic moment, 
e = mv 2 /2 is particle energy, Fq is a Maxwellian, ip is the 
electrostatic potential, Bq is the magnitude of the back- 
ground magnetic field, T is the background temperature, 
q is particle charge, and f2o = qBa/mc. The gyrokinetic 
equation governing the evolution of h is given by 
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where b = B0/-B0, is the drift velocity of guiding 
centers, x = — v-A/c, A is the vector potential, {a, b} 



is the Poisson bracket of a and b, 
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is the gyroaver- 



age of a at constant R, and (C[h]) R is the gyroaveraged 
collision operator. 

For (C[h]) R , we restrict our attention to the model 
collision operator presented in Paper I. We work within 
the framework of the continuum gyrokinetic code GS2, 5 
which assumes periodicity in the spatial directions per- 
pendicular to Bo in order to reduce the simulation vol- 
ume to a thin flux tube encompassing a single magnetic 
field line. Consequently, we require a spectral represen- 
tation of (C[h]) R : 



(C[h]) n = J2^ n C G K[h k }, 



(5) 



where k is the perpendicular wave vector. For conve- 
nience, we reproduce the expression for the same-species 
part of CcKihu] from Paper I in operator form: 

C GK [h k ) = L[h k ]+D[h k ]+U L [h k ] + U D [h k ]+E[h k ], (6) 

where 
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are the gyroaveraged Lorentz and energy diffusion oper- 
ators (which together form the test-particle piece of the 
linearized Landau operator, as shown in Refs 1291 and |3TJ]| . 
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are the gyroaveraged momentum-conserving corrections 
to the Lorentz and energy diffusion operators, and 



E[h k ] = v E v 2 J {a)F a 
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is the gyroaveraged energy-conserving correction (the 
conserving terms are an approximation to the field- 
particle piece of the linearized Landau operator). The 
electron collision operator has the following additional 
term to account for electron-ion collisions: 
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where £ = v\\/v is the pitch angle, a = kvj_/flQ, J and 
Ji are Bessel functions of the first kind, v t h = y/2To/m is 
the thermal velocity, and My [/i£,k] is the perturbed parallel 
ion flow velocity Expressions for the velocity-dependent 
collision frequencies vd, Av, i/y, and ve are given in Pa- 
per 1 (which follows the notation of Ref. l3"Tj) . 

Having specified the form of our collision operator, we 
now discuss some of its fundamental properties that guide 
our choice of numerical implementation. 



A. Collision operator amplitude 

Even when the collisionality approaches zero, C G K[h k ] 
can have appreciable amplitude. There are two reasons 
for this: first, the velocity dependence of vp, ve, and 
Av is such that each go to infinity as v — > (so low- 
velocity particles arc always collisional) ; and second, we 
expect the distribution function to develop increasingly 
smaller scales in v and £ as collisionality decreases, so 
that the amplitude of the terms proportional to d 2 h/d(, 2 
and d 2 h/dv 2 may rema in approxi mately constant!^! (i.e. 
C GK [h k } 7^ as v -> 0) |iu | ii | i2 | i4 | Tne fact that C GK [h k ] 



can be quite large even at very low collisionalities means 
that it should be treated implicitly if one wants to avoid a 
stability limit on the size of the time step, At. In Section 
III, we describe our fully implicit implementation of the 
collision operator. 



B. Local moment conservation 

Since collisions locally conserve particle density, mo- 
mentum, and energy, one would like these properties to 
be guaranteed by the discrete version of the collision op- 
erator. Mathematically, this means that the density, 
momentum, and energy moments of the original (un- 
gyroaveraged) collision operator must vanish (for same- 
species collisions). However, the non-local nature of the 
gyroaveraging operation introduces finite Larmor radius 
(FLR) effects that lead to nonzero values for the analo- 
gous moments of (C[/i]) R . Since this is the quantity we 
employ in gyrokinetics, we need to find the pertinent re- 
lations its moments must satisfy in order to guarantee 
local conservation properties. 

This is accomplished by Taylor expanding the Bessel 
functions Jo and J\. In particular, one can show thaiPsl 

1 
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where (.) r denotes a gyroaverage at fixed r, C GK [h k ] is 
the operator of Eqn. (|6| with kp = (neglecting FLR 
terms but retaining nonzero subscripts k for h), and Tq 
is the collisional flux of number, momentum, and energy 
arising from FLR terms. Consequently, the density, mo- 
mentum, and energy moments of the gyrokinetic equa- 
tion can be written in the conservative form 
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where M. = (Sn nSu^ Sp) represents the perturbed 
number, momentum, and energy densities, the super- 
script T denotes the transpose, and Tm contains both 
the collisional flux, Tq, and the flux arising from all other 
terms in the gyrokinetic equation (for a more detailed 
discussion, see Paper I). Thus, local conservation prop- 
erties are assured in gyrokinetics as long as the density, 
momentum, and energy moments of C GK [h k ] vanish: 
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(15) 



We describe how this is accomplished numerically in Sec- 
tion III. 
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C. //-Theorem 

In contrast with local conservation properties, the 
statement of the //-Theorem is unmodified by gyroav- 
eraging the collision operator. Defining the entropy as 
5 = —/In/, Boltzmann's //-Theorem tells us 



35 
dt 



V 



d 3 v ln[f}C[f] > 0, (16) 



where V = J d 3 r and the double integration spans phase 
space (the velocity integration is taken at constant par- 
ticle position r). Expanding the distribution function as 
before, we find to lowest order in the gyrokinetic ordering 
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d 3 v —C[h] < 0. 
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Changing variables from particle position r to guiding 
center position R, we obtain 



d 3 v 



d 3 R h 
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(C[k]) R < 0; 



(18) 



where now the velocity integration is taken at constant 
R. In this case, the non-locality of the gyroaveraging 
operation leads to no modification of the //-Theorem be- 
cause of the definition of entropy as a phase-space aver- 
aged quantity (as opposed to local conservation proper- 
ties, which involve only velocity-space averages). There- 
fore, one can easily diagnose entropy generation and test 
numerical satisfaction of the //-Theorem in gyrokinetic 
simulations, as we show in Sec. IV. 



III. NUMERICAL IMPLEMENTATION 

It is convenient for numerical purposes to separately 
treat collisional and collisionlcss physics. Thus, we begin 
by writing the gyrokinetic equation in the form 



dt 



= C GK [h k ] +A[h k ], 



(19) 



where -A[/ik] represents the rate of change of h k due to 
the collisionlcss physics. In order to separate these terms, 
we utilize Godunov dimensional splitting,^ which is ac- 
curate to first order in the timestep At: 



At 

K +1 - h 
At 



(20) 
(21) 



where n and n + 1 are indices representing the current 
and future time steps, and h k is defined by Eqn. ( 20 1 
it is the result of advancing the collisionless part ofthc 
gyrokinetic equation. With h ^. th us given, we restrict 
our attention to solving Eqn. pTl. For notational con- 
venience, we suppress all further k subscripts, as we will 
be working exclusively in fc-space. 



As argued in Section II, we must treat the collision 
operator implicitly to avoid a stability limit on the size 
of At. We use a first order accurate backward-difference 
scheme in time instead of a second order scheme (such 
as Crank-NicholsorPSf) because it is well known that 
the Crank-Nicholson scheme introduces spurious behav- 
ior in solutions to diffusion equations when taking large 
timesteps (and because Godunov splitting is only first 
order accurate for multiple splittings, which will be in- 
troduced shortly). 

With this choice, h n+1 is given by 



(1-AtCWr 1 ^*. 



(22) 



In general, Cgk is a dense matrix, with both energy and 
pitch-angle indices. Inversion of such a matrix, which 
is necessary to solve for h n+1 in our implicit scheme, is 
computationally expensive. We avoid this by taking two 
additional simplifying steps. First we employ another 
application of the Godunov splitting technique, which, 
combined with the choice of a (£, v) grid in GS2J31 al- 
lows us to consider energy and pitch-angle dependence 
separately!^ 

h ** = [l_At(L + U L )]- 1 h* (23) 
h n+1 = [1- At(D + U D + E)}' 1 h**, (24) 

The h* and h** are vectors whose components are the 
values of h at each of the (£, v) grid points. In Eqn. (23 1 
we order the components so that 



h — (hii,h 2 i, fiNi, h\2, HnmY 



(25) 



where the first index represents pitch-angle, the second 
represents energy, and N and M are the number of pitch- 
angle and energy grid points, respectively. This allows for 
a compact representation in pitch-angle. When solving 
Eqn. (24), we reorder the components of h so that 



h = (hi 1 ,hi2,—,hiN,h 2 i,— ,/ijvm) , (26) 
allowing for a compact representation in energy. 



A. Conserving terms 

The matrices 1 — AtL and 1 — AtD are chosen to be 
tridiagonal by employing three-point stencils for finite 
differencing in £ and v. This permits computationally 
inexpensive matrix inversion. However, the full matri- 
ces to be inverted include the momentum- and energy- 
conserving operators, U and E, which are dense matrices. 
We avoid direct inversion of these matrices by employing 
the Sherman-Morrison formulaf 34 * 35 ! which gives x in the 
matrix equation A/x = b, as long as M can be written 
in the following form: 



M = A 



(27) 
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where £g> is the tensor product . The solution is then given 

by 



V • y 



1 



(28) 



where y = A _1 b, z = A _1 u, and the dot products repre- 
sent integrals over velocity space. If A^ 1 is known or eas- 
ily obtainable (as in our case), this formulation provides 
significant computational savings over the straighforward 
method of directly inverting the dense matrix M. 

Details of the application of the Sherman- Morrison for- 
mula to Eqns. (23 1 and (24) are given in Appendix A. 
Here, we state the main points. The matrix operators 
L and D are to be identified with A, and the integral 
conserving terms U and E can be written in the form of 
the tensor product, u£g)v. Identifying h** and h n+1 with 
x, we find that multiple applications of the Sherman- 
Morrison formula give 
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(29) 

w , (30) 
w . (31) 



The quantities Vq, Vi, V2, Zq, So, wo, and yo are specified 
in Table[l]in Appendix A. With the exception of yo, each 
of these quantities is time-independent, so they need be 
computed only once at the beginning of each simulation. 
Consequently, inclusion of the conserving terms in our 
implicit scheme comes at little additional expense. 



We note that when Eqn. (29) is applied to comput 



ing the inverse matrix in Eqn. (23), the corresponding 



v 2 is nonzero only for the electron collision operator. 
This term arises by using the parallel component of Am- 
pere's law to rewrite the electron-ion collison operator of 
Eqn. (Il2|) as 
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where e is the magnitude of the electron charge, u» [h e ] 
is the parallel component of the electron fluid velocity, 
and no, e is the equilibrium electron density. For electron 
collisions, the Ay term is absorbed into h* so that we use 
the modified quantity 



v%At 



ck 2 v\\ 



2nev 2 h no, 



-Ay J (a e )F 0t , 



(33) 



when applying the Sherman-Morrison formula, where Ay 
from the n+1 time level is used (for details on the implicit 
calculation of Ay, see Ref. |5J. 

B. Discretization in energy and pitch angle 

We still must specify our choice of discretization for 
Cqk- Ideally, we would like the discrete scheme to guar- 
antee the conservation properties and ff-Theorem asso- 
ciated with C. As discussed in Sec. II, the former is 
equivalent to requiring t hat the kp — component of 
Cgk, Cq K , satisfy Eqn. (151. We now proceed to show 
that this requirement is satisfied by carefully discretizing 
the conserving terms and by employing a novel finite dif- 
ference scheme that incorporates the weights associated 
with our numerical integration scheme. 

We begin by writing C GK [h] for same-species collisions: 
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J d v i>DV\\h 
J d 3 v v D vjF 



(34) 



With C GK thus specified, we now consider numerical 
evaluation of the relevant moments of Eqn. (34). To sat- 
isfy number conservation (J d 3 v Cq K [Ii] =0), velocity 
space integrals of each of the terms in Eqn. (34 1 should 



vanish individually. For integrals of the first two terms 
to vanish, we require a finite difference scheme that sat- 
isfies a discrete analog of the Fundamental Theorem of 
Calculus (i.e. conservative differencing); for the last three 
terms, we must have a discrete integration scheme satisfy- 



ing J d 3 v ise>v\\F = J d 3 v Avv\\Fo = J d 3 v uev 2 F q = 0. 
The requirement that J d 3 v ^.Dfy-fo = / d 3 v AvvhFq = 
is satisfied by any integration scheme with velocity space 
grid points and associated integration weights symmetric 
about uy = 0, which is true for the (£, v) grid used in GS2. 
By substituting for ue everywhere using the identity 
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\v 5 F ) , 



(35) 
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the other integral constraint (J d 3 v vev 2 Fq = 0) reduces 
to the requirement that finite difference schemes must 
satisfy the Fundamental Theorem of Calculus. 

Parallel momentum conservation (J d 3 v v\\Cq K [}i] = 
0) introduces the additional requirements that: 



d 3 v 



and 
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d 3 v Avv 2 F 



fd 3 vAvv\\h (37) 
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If the finite difference scheme used for all differentiation 
possesses a discrete version of integration by parts (upon 
double application), then Eqns. (36 1 and (37 1 are numer- 
ically satisfied as long as: Vnvoh in the second term of 
Eqn. d36|) is expressed in the form 



v\\V£>h 



1 



d_ 



dvn 
~dj 



v D h, 



(38) 



Av on the righthand side of Eqn. (37 1 is expressed using 
the identity 



2Avv 3 Fn 



— (v v A F — 
dv V ° dv 



(39) 



and all integrals are computed using the same numerical 
integration scheme (if analytic results for the integral de- 
nominators in terms 3 and 4 of Eqn. (34 1 are used, then 



the necessary exact cancellation in Eqns. (36) and (37) 
will not occur). 

The only additional constraint imposed by energy con- 
servation (J d 3 v v 2 CQ K [h] = 0) is that the form of 
Eqn. p5|) be slightly modified so that 



u E v 2 F 
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dv 2 
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(40) 



which still satisfies the number conservation contraint. 



Using the forms given by Eqns. ( 38 1-( 40 ) , conservation 



properties are guaranteed as long as one employs a finite 
difference scheme for pitch-angle scattering and energy 
diffusion that satisfies discrete versions of the Fundamen- 
tal Theorem of Calculus and integration by parts. 

For the case of equally spaced grid points in v and 
£, there is a straightforward difference scheme, accurate 
to secord order in the grid spacing, that satisfies both 



.36 



requirements 

d ^dh ^ Gj+1/2 (hj+i - hj) - Gj-i/2 (hj - foj-i) 
dx dx Ax 2 

(41) 

where X is a dummy variable representing either v or £, 
Ax is the grid spacing, hj is the value of h evaluated at 



FIG. 1: (Left): Solid line indicates the scaling of the leading 
order error, averaged over all grid points, of the conservative 
finite difference scheme for a Gauss-Legendre grid (the grid 
used in GS2). The slope of the dotted line corresponds to a 
first order scheme. (Right): factor by which the conservative 
finite difference scheme of Eqn. ( 42 ) amplifies the true collision 



operator amplitude at the boundaries of the Gauss-Legendre 
grid. 



the grid point Xj, £j±i/2 = ( x j + x j±i)/%> an d G is either 
1 — £ 2 (for pitch-angle scattering) or i/nv Fq (for energy 
diffusion). However, in order to achieve higer order ac- 
curacy in the calculation of the velocity space integrals 
necessary to obtain electromagnetic fields, GS^S! and a 
number of other gyrokinetic codeipZl use grids with un- 
equal spacing in v and £ and integration weights that are 
not equal to the grid spacings. 

Given the constraints of a three-point stencil on an un- 
equally spaced grid, we are forced to choose between a 
higher order scheme (a second order accurate scheme can 
be obtained with compact differencing tSSl as described in 
Appendix B) that does not satisfy our two requirements 
and a lower order scheme that does. Since our analytic 
expression for Cqk wa s designed in large part to sat- 
isfy conservation properties (and because the conserving 
terms are only a zeroth order accurate approximation to 
the field-particle piece of the linearized Landau opera- 
tor^), we choose the lower order scheme, given here, as 
the default: 

9 r dh ~ 1 (r Z h J r h. Z h i~ l 

(j-j + i/2 Lr,-_i/2 

dx dx Wj \ x j+i ~ x j x j ~ x j-i 

where Wj is the integration weight associated with Xj. 

Defining * = Gh! , with the prime denoting differen- 
tiation with respect to x, Taylor series can be used to 
show 



ttj+i/2 - ttj-1/2 = Axj | Q ( (Ax); 



where Ax,- 



C J + l/2 



Cj-l/2- 



(43) 



With the exception of 



pitch angles corresponding to trapped particles J*^ 1 the 
grid points {xj} and associated integration weights {vjj} 
in GS2 are chosen according to Gauss-Legendre quadra- 
ture rules For this case, we show numerically in Fig. [l] 
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FIG. 2: Plots showing evolution of the perturbed local density, parallel momentum, and energy over fifty collision times. 
Without the conserving terms (TsJ)— < 11 ) , both parallel momentum and energy decay significantly over a few collision times 
(long dashed lines). Inclusion of conserving terms with the conservative scheme detailed in Sec. Ill leads to exact moment 
conservation (solid lines). Use of a non-conservative scheme leads to inexact conservation that depends on grid spacing (short 
dashed lines). 



that 



N 



and 



1 N 

-T 



= 1+0 



max 

j=2,...,N-l 



= 1+0 



1 A Xj 








Wj 





(44) 



(45) 



where N is the number of grid points in x. 

The boundary points (J = 1, N) are excluded from the 
max operator above. This is because Ax/w (the factor 
multiplying ^ in Eqn. (43 1 ) converges to approximately 
1.2 for the boundary points as the grid spacing is de- 
creased (Fig. [I]). For the energy diffusion operator, we 
can make use of the property that G(x) = G(x)' = at 
x = and x = oo to show 



± 



i±i/2 




(46) 



with the plus sign corresponding to 7 = 1 and the minus 
sign to j = N. This is not true for the Lorentz operator, 
so we are forced to accept an approximately twenty per- 
cent magnification of the Lorentz operator amplitude at 
£ = ±1 and at the trapped-passing boundaries. We find 
that this relatively small error at the boundaries has a 
negligible effect on measureable (velocity space averaged) 
quantities in our simulations. 



IV. NUMERICAL TESTS 

We now proceed to demonstrate the validity of our col- 
lision operator implementation. In particular, we demon- 



strate conservation properties, satisfaction of Boltz- 
mann's i?-Theorem, efficient smoothing in velocity space, 
and recovery of theoretically expected results in both col- 
lisional (fluid) and collisionlcss limits. While we do not 
claim that the suite of tests we have performed is exhaus- 
tive, it constitutes a convenient set of numerical bench- 
marks that can be used for validating collision operators 
in gyrokinetics. 



A. Homogeneous plasma slab 

We first consider the long wavelength limit of a homo- 
geneous plasma slab with Boltzmann electrons and no 
variation along the magnetic field (fcy — 0). The gyroki- 
netic equation for this system simplifies to 



d(Sf) 
dt 



C° GK [h], 



(47) 



which means local density, momentum, and energy 
should be conserved. In Fig. [2] we show numerical results 
for the time evolution of the local density, momentum, 
and energy for this system. 

Without inclusion of the conserving terms j9|-( 11 1, we 
see that density is conserved, as guaranteed by the con- 
servative differencing scheme, while the momentum and 
energy decay away over several collision times. Inclusion 
of the conserving terms provides us with exact (up to nu- 
merical precision) conservation of number, momentum, 
and energy. To illustrate the utility of our conservative 
implementation, we also present results from a numer- 
ical scheme that does not make use of Eqns. ( 38 1-( 40 > 



and that employs a finite difference scheme that does not 
possess discrete versions of the Fundamental Theorem of 
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FIG. 3: Plot of the evolution of entropy generation for the 
homogeneous plasma slab over twenty collision times. Our 
initial distribution in velocity space is random noise, and we 
use a grid with 16 pitch angles and 8 energies. The entropy 
generation rate is always nonnegative and approaches zero in 
the long-time limit. 



FIG. 4: Evolution of |Jm| for the electromagnetic plasma slab 
with (3 = 10 -4 , k y pi = 0.1, and u e t = 10k^v t h,i- Inclusion of 
the ion drag term in the electron-ion collision operator leads 
to the theoretically predicted damping rate for the parallel 
current given in Eqn. ( 49 1 . Without the ion drag term, the 
parallel current decays past zero (at t ~ 22) and converges to 
a negative value as the electron flow damps to zero. 



Calculus and integration by parts. Specifically, we con- 
sider a first order accurate finite difference scheme similar 
to that given by Eqn. ( |42| , with the only difference be- 
ing that the weights in the denominator are replaced with 
the local grid spacings. In this case, we see that density, 
momentum, and energy are not exactly conserved (how 
well they are conserved depends on velocity space reso- 
lution, which is 16 pitch angles and 16 energies for the 
run considered here). 

The rate at which our collision operator generates en- 
tropy in the homogenous plasma slab is shown in Fig. [3] 
As required by the iJ-Thcorem, the rate of entropy pro- 
duction is always nonnegative and approaches zero in the 
long-time limit as the distribution function approaches a 
shifted Maxwellian. We find this to hold independent 
of both the grid spacing in velocity space and the ini- 
tial condition for the distribution function (in Fig. [3j the 
values of /i(£, v) were drawn randomly from the uniform 
distribution on the interval [—1/2, 1/2]). 



B. Resistive damping 



gleet d(5f)/dt. However, since An ~ k~ 2 , and we are 
considering t « 1, dA^/dt must be retained. The re- 
sulting electron equation is of the form of the classical 
Spitzer problem (see, e.g., Ref. 40): 



^GK 



[h c 



T e c dt 



(48) 



The parallel current evolution for this system is given by 
JM) = J,|(t = 0)e-* (49) 



1.98r e ri e e 2 /m e c 2 
dy/w/Avei is the 



where r\ = l/c|| is the resistivity, o\\ 
is the Spitzer conductivity, and r e 
electron collision time. 

We demonstrate that the numerical implementation of 
our operator correctly captures this resistive damping in 
Figs. [4] and [5] We also see in these figures that in the 
absence of the ion drag term from Eqn. (l2"| , the electron 
flow is incorrectly damped to zero (instead of to the ion 
flow), leading to a steady-state current. 



We now modify the system above by adding a finite A\\ . 
From fluid theory, we know that collisional friction be- 
tween electrons and ions provides resistivity which leads 
to the decay of current profiles. Because the resistive 
time is long compared to the collision time, one can ne- 



C. Slow mode damping 

We next consider the damping of the slow mode in a 
homogenous plasma slab as a function of collisionality. 
In the low k±pi, high [3i limit, one can obtain analytic 
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FIG. 5: Evolution of perturbed parallel flow for the elec- 
tromagnetic plasma slab with (3 = 10 -4 , k y pi — 0.1, and 
v e i — lOfcyUti. Without inclusion of the ion drag term in 
Eqn. (121, the electron flow is erroneously damped to zero 



(instead of to the ion flow). 



expressions for the damping rate in both the collisional 
(k»\ m fp <C 1) and collisionlcss (k^\ m f p ^> 1) regimes, 
where X m fp is the ion mean free path (see e.g. Ref. [T2]) . 
The expressions are 



±knv A \ I 



for fc||A m / p <C 1, and 



2v A 



\k\\\v A 



(50) 



(51) 



FIG. 6: Damping rate of the slow mode for a range of col- 
lisionalities spanning the collisionless to strongly collisional 
regimes. Dashed lines correspond to the theoretical predic- 
tion for the damping rate in the collisional (fc||A m / p <C 1) and 
collisionless (fc||A m / p 2> 1) limits. The solid line is the result 
obtained numerically with GS2. Vertical dotted lines denote 
approximate regions (collisional and collisionless) for which 
the analytic theory is valid. 



behavior in the k\\\ m f p <C 1 limit (damping rate pro- 
portional to the correct collisional damping in the 
k\\\ m f p ~ 1 limit (damping rate inversely proportial to 
and the correct collisionless (i.e. Barnes) damping 
in the k\\\ m f p 1 limit. 



D. Electrostatic turbulence 



for knX m f p 1. Here, v A = v t h,i/\fWi is the Alfven 
speed, and i/y^ is the parallel ion viscosity, which is in- 
versely proportial to the ion-ion collision frequency, v>u: 
jva. As one would expect, the damping in 



J th,i 



the strongly collisional regime [Eqn. (50)] is due primar- 



ily to viscosity, while the collisionless regime [Eqn. (51 1] 
is dominated by Barnes damping.^ 

In Fig. [6] we plot the collisional dependence of the 
damping rate of the slow mode obtained numerically us- 
ing the new collision operator implementation in GS2. In 
order to isolate the slow mode in these simulations, we 
took ip = j4|| = Sn e = and measured the damping rate 
of SBu . This is possible because 8B\\ effectively decouples 
from ip and A\\ for our system, and Sn e can be neglected 
because 3> lP^ We find quantitative agreement with 
the analytic expressions (50) and (51) in the appropri- 



ate regimes. In particular, we recover the correct viscous 



Finally, we illustrate the utility of our collision opera- 
tor in a nonlinear simulation of electrostatic turbulence in 
a Z-pinch field configuration.^ We consider the Z-pinch 
because it contains much of the physics of toroidal con- 
figurations (i.e. curvature) without some of the complex- 
ity (no particle trapping). At relatively weak pressure 
gradients, the dominant gyrokin etic linear ins tability in 
the Z-pinch is the entropy mode,' 43 ' 44 ' 45 ' 46 ' 47 ' 48 ' which is 
nonlincarly unstable to secondary instabilities such as 
Kelvin- Helmholtz. 49 

In previous numerical investigations of linear^Sl and 
nonlinear^ plasma dynamics in a Z-pinch, collisions were 
found to play an important role in the damping of zonal 
flows and in providing an effective energy cutoff at short 
wavelengths. However, as pointed out in Ref. |49j the 
Lorentz collision operator used in those investigations 
provided insufficient damping of short wavelength struc- 
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FIG. 7: Evolution of ion particle and heat fluxes for an elec- 
trostatic, 2-species Z-pinch simulation. We are considering 
R/ L„ — 2.0 and va — 0.01vth,i/R- The particle flux is indi- 
cated by the solid line and is given in units of (p/R)no,iVth,i- 
The heat flux is indicated by the dashed line and is given in 
units of (p/7?)no,i^,i- We see that a steady-state is achieved 
for both fluxes without artificial dissipation. 



tures to obtain steady-state fluxes. Consequently, a 
model hyper- viscosity had to be employed. 

We have reproduced a simulation from Ref. HH1 us- 
ing our new collison operator, and we find that hyper- 
viscosity is no longer necessary to obtain steady-state 
fluxes (Fig. [7}. This can be understood by examining 
the linear growth rate spectrum of Fig. [8] We see that 
in this system energy diffusion is much more efficient at 
suppressing short wavelength structures than pitch-angle 
scattering. Consequently, no artificial dissipation of short 
wavelength structures is necessary. 



V. SUMMARY 

In Sec. I we proposed a set of key properties that an 
ideal dissipation scheme for gyrokinetics should satisfy. 
Namely, the scheme should: limit the scale size of struc- 
tures in phase space in order to guarantee the validity of 
the gyrokinetic ordering and to provide numerical res- 
olution at reasonable expense; conserve particle num- 
ber, momentum, and energy; and satisfy Boltzmann's 
ff-Theorem. While commonly employed simplified col- 
lision operators or hyperviscosity operators may be ade- 
quate for some calculations,^ it is important to be able 
to use the more complete collision operator described in 
this paper, which preserves all of these desirable dissipa- 



FIG. 8: Linear growth rate spectrum of the entropy mode in 
a Z-pinch for R/L n = 2.0, where R is major radius and L n 
is density gradient scale length. The solid line is the colli- 
sionless result, and the two dashed lines represent the result 
of including collisions. The short dashed line corresponds to 
using only the Lorentz operator, while the long dashed line 
corresponds to using our full model collision operator. Both 
collisional cases were carried out with va = Q.Qlv t h,i/ R- 



tion properties. 

In Sec. II we presented the model collision operator 
derived in Paper I and discussed some of its features 
that strongly influence our choice of numerical imple- 
mentation. In particular, we noted that local conserva- 
tion properties are guaranteed as long as the (1, v», v 2 ) 
moments of the kp — component of the gyroaveraged 
collision operator vanish. Further, we argued that the 
collision operator should be treated implicitly because in 
some regions of phase space, its amplitude can be large 
even at very small collisionalities. 

Our numerical implementation of the collision oper- 
ator was described in Sec. III. We separate collisional 
and collisionlcss physics through the use of Godunov di- 
mensional splitting and advance the collision operator 
in time using a backwards Euler scheme. The test par- 
ticle part of the collision operator is differenced using a 
scheme that possesses discrete versions of the Fundamen- 
tal Theorem of Calculus and integration by parts (upon 
double application). These properties are necessary in 
order to exactly satisfy the desired conservation prop- 
erties in the long wavelength limit. The field particle 
response is treated implicitly with little additional com- 
putational expense by employing repeated application of 
the Sherman-Morrison formula, as detailed in Appendix 
A. 
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In Sec. IV we presented numerical tests to demon- 
strate that our implemented collision operator possesses 
the properties required for a good gyrokinetic dissipation 
scheme. In addition to these basic properties, we showed 
that the implemented collision operator allows us to cor- 
rectly capture physics phenomena ranging from the colli- 
sionless to the strongly collisional regimes. In particular, 
we provided examples for which we are able to obtain 
quantitatively correct results for collisionless (Landau or 
Barnes), resistive, and viscous damping. 

In conclusion, we note that resolution of the collision- 
less (and collisional) physics in our simulations was ob- 
tained solely with physical collisions; no recourse to any 
form of artificial numerical dissipation was necessary. 
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VI. APPENDIX A: SHERMAN-MORRISON 
FORMULATION 

The repeated application of the Sherman-Morrison for- 
mula considered here is an extension of the scheme pre- 
sented in Tatsuno and Dorland.^ Throughout this cal- 
culation, we adopt general notation applicable to both 
Eqns. ( 23 1 and p4| and provide specific variable defini- 
tions in Table [I| Both Eqns. ([23) and ((24 1 can be written 
in the form 



Ax 



(52) 



Because (14) and (£) are integral operators, we can write 
them as tensor products so that 



Uq ® V + Ui ® Vi + u 2 <£> v 2 . 



A= A 
We now define 

Ai = A + u <g) v , A 2 = Ax + ui <g> vi, 
so that 

(A 2 + u 2 <g)v 2 )x = b. 



(53) 



(54) 



(55) 



Applying the Sherman-Morrison formula to this equa- 
tion, we find 



y2 



V2 ■ Y2 

1 + v 2 • z 2 



Z2. 



(56) 



variable 


L 


D 


A 


l-At(L + U L ) 


1 — At (D + Ud + E) 


X 


h** 


h n+1 


b 


h* 


h** 


A 


1 - AtL 


1 - AtD 


Vo 


v D v±J 1 


-AwiJi 


Vl 


V D V\\J 


— Avvu Jo 


v 2 


H|| (electrons) 
(ions) 


f E v 2 J 


U() 


— Atv /d u 


-Atvo/d u 


Ul 


-Atvi/du 


—Atvi/du 


u 2 


— AtVpV\\/d q (electrons) 
(ions) 


-Atv 2 /d q 


d u 


/ d 3 v vdVuFo 


jd 3 v Ai/ujjFo 


d q 




/ d 3 v v E v i F 



TABLE I: Sherman-Morrison variable definitions for Lorentz 
and energy diffusion operator equations 



where A 2 y 2 = b and A 2 z 2 = u 2 . 

Applying the Sherman-Morrison formula to each of 
these equations gives 



Y2 = yi - 
z 2 = zi - 



vi - yi 

1 + V! • Wj 
Vl • Zl 



Wi 



1 + Vl • Wi 



Wl, 



(57) 
(58) 



where Aiyi = b, AiWi = Ui, and AiZi = u 2 . A final 
application of Sherman- Morrison to these three equations 
yields 



v • yo 

yi = yo - \ . s 



Wi = Wo 

zi = z - 



1 + v • s 
vp ■ w 

1 + v • s 
v ■ z 



S() 



1 + v • s 



So- 



(59) 
(60) 
(61) 



where A y = b, A s = u , A w = Ui, and A z = 
u 2 . 

We can simplify our expressions by noting that Vo,i,2 
and Uo,i,2 have definite parity in v\\. A number of inner 
products then vanish by symmetry, leaving the general 
expressions 



y2 = yo 



z 2 = z 



vq • yo 
1 + v • s 

vq ■ z 
1 + v • s 



s - 
so - 



vi • yo 

1 + Vl • w 

vi ■ z 
1 + vi • w 



w (62) 
w Q . (63) 
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VII. APPENDIX B: COMPACT 
DIFFERENCING OF THE TEST-PARTICLE 
OPERATOR 

In this appendix, we derive a second order accurate 
compact difference scheme for pitch-angle scattering and 
energy diffusion on an unequally spaced grid. The higher 
order of accuracy of this scheme is desirable, but it does 
not possess discrete versions of the Fundamental The- 
orem of Calculus and integration by parts when used 
with Gauss-Legendre quadrature (or any other integra- 
tion scheme with grid spacings unequal to integration 
weights). Consequently, one should utilize this scheme 
only if integration weights and grid spacings are equal, 
or if exact satisfaction of conservation properties is not 
considered important. 

For convenience, we begin by noting that Eqns. (23) 
and ( 24 I can both be written in the general form 



dh\ 

at) ( 



= H (Gti)' + S = H(G'ti + Gh") + S, (64) 



where: for the Lorentz operator equation (23) we iden- 
tify H = 1, G = v D {l-e)/2, S = U L [h] - 



k v Vd (1 + £ ) /j/4f2o ; an d the prime denotes differ- 
entiation with respect tof; and for the energy diffu- 
sion operator equation (24), we identify H = l/2v 2 F 
G = ^v 4 F Q ,S 



U D [h] + E[h] - k 2 v 2 v\\ (1 - ( 2 ) h/iQl 
and the prime denotes differentiation with respect to v. 
Here, the h we are using is actually normalized by Fq. 

Employing Taylor Series expansions of h, we obtain 
the expressions 



ti 



S 2 . (fej+i - hj) + 5+ (hj - fej-i) 
5+5- (5+ + 6-) 



+ 0[5\ (65) 



and 



ti' = 2 



5- (h i+ i - hj) - 5+ (hj - hj-i) 5- - 5 + 
6+6- (5+ + 5_) 3 



h'l'+0[5 2 
(66) 

where i denotes evaluation at the velocity space gridpoint 
Xi, and 5± = \xi±\ — Xi\ (here i is a dummy variable 
representing either £ or v). In order for the h'( expression 
to be second order accurate, we must obtain a first order 
accurate expression for h"' in terms of hi, h[, and h". We 



accomplish this by differentiating Eqn. (64) with respect 
to x: 



dti 
~dt 



= H' (Gti)' + H (Gti)" + S' (67) 



and rearranging terms to find 



h? = 



1 



HiG 
Hi (G'!h 



dh,, 
~dt 



- Hi (G'jh'i + Gih'l) 

c 

2G' l til)-S'] +0[5], 



(68) 



where, unless denoted otherwise, all quantities are taken 



at the n + 1 time level. Plugging this result into Eqn. 66 



and grouping terms, we have 

5- -5. 



mK = 



">* I ~ ti l G i - H l G i I — bi 



iHiGi 

5- (h l+ i - ti) - 5+ (hj - hj-j) 
6+6- (5+ + 5-) 



0[5 2 



(69) 



where ^ = 1 + (5+ - 5-) (H[Gi + 2H i G' i ) /3HiG h and 



we have taken (dh/dt)c 



(ti' 



h™)/At. Using 



Eqn. (67) and the above result in Eqn. (|64|), we find 

5+ + 5. 



dh 



c 



= ti, [ HG 1 , - 



3^ 4 



1 
At 



tji r^f I I fill 
— ±1 % \JJ . ( — STL I Kjr ^ 



HjGj ( 6- (hj+i - ti) - 5+ (hj ■ 
Hi \ 6+5-(6+ + 5_) 



l) 



J + 



3HGi 



(KY 

At 



s' 



S q + 0[5 2 



(70) 



This is the general compact differenced form to be used 
when solving Eqns. (23 1 and (24 1. 



In order to illustrate how compact differencing affects 
the implicit solution using Sherman-Morrison, we present 
the result of using the particular form of S for energy 
diffusion in Eqn. (70): 



h? 



H+l 



At 



6+ (6+ 
hi-i 



2H t G, 

5-) { m 

( 2HiGi 



5- (5+ + S-) V H 
hi ( ^H^G^ 



+ 5-0 




+ (5+-5-)Q l -5+5_Kv t 



- W 



6-) 



5- (6+ + 5-) 



+ U ll V ll +U+V+ + U q q + 0[5 2 }, 

(71) 



where 



(Ti = (6+ - 5-) /3m 

k = k 2 v 2 th (i - e) m 2 

Q = HiG't - a, (1/At - H'fi'i ~ HG'l + Kv s ) 
v s = i>\\v 2 th /2v 2 
A = A + a A! 

^-L,||,g = u 0,l,2 
KL,||,g = v 0,l,2 

with u and v given in Table [TJ 

We see that the only significant effects of compact dif- 
ferencing on numerical implementation are: modification 
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of h** in Eqn. (24) to reflect the h n terms on the right- 
hand side of Eqn. (71 1; and modification of the Un , U±, 



and U q terms that appear in Sherman-Morrison (uo,i,2 
from Appendix A) to include an additional aU' term. 
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